#install libraries
library(ape)
library(stats)
library(geiger)
library(ade4)
library(apTreeshape)
library(methods)
#library(phybase)
library(phylobase)




#read table with data:
mydata=read.table("allCOdata.txt", sep="\t", header=TRUE)
##put seq_id as row.names
#names_mydata=gsub(">","",mydata$seq_names)
row.names(mydata)=mydata$seq_names


#read tree Newick format (branch lenghts and nodes identified)
mytree=read.tree("CO1_20100809_FastTr")


#plot tree
pdf('mytree_plot_CO1.pdf')
plot.phylo(mytree)
dev.off()


#check taxa in table are the same as in the tree
name.check(mytree,mydata)


#calculate patristic distances matrix from the tree (sum of the lengths of the branches that link two nodes in a tree)
myPatrisDis=cophenetic(mytree)
write.matrix(myPatrisDis,"CO1PatrisDist.txt")


#clustering hierarchically (threshold t% for clusters defined by user)from the patristic distance matrix
##hclust() needs an object of class 'dist'
myMatrix=as.dist(myPatrisDis)

## single adopts a 'friends of friends'strategy, ie, the distance between 2 clusters is the minimum of the distances between elements of the 2 clusters
clusterT=hclust(myMatrix, method="single")

##save plot as pdf in the present working directory
pdf('clusterT_plot_CO1.pdf')
plot(clusterT)
dev.off()

##clusters are formed using the threshold value=h (cutree needs an hclust class object)
thr=c(0.01,0.02,0.05,0.07,0.1)
print(thr)
threshold=list()
for (x in 1:length(thr))
{
threshold[[x]]=cutree(clusterT, h=thr[x])
names(threshold)[x]=paste("threshold_",thr[x],sep="")
}
print(threshold)
threshold.df=data.frame(threshold)

##to print boxes in the dendrogram around the clusters obtained...it's not working!!
#rect.hclust(clusterT, h=0.01, border="red")




#import to mydata the threshold column (with the different clusters obtained at threshold t%)

##adding a 'seq_names' column with the row.names to threshold.df
new_thr=transform(threshold.df,seq_names=row.names(threshold.df))
new_thr$seq_names=as.character(new_thr$seq_names)   #maybe not needed!

##add the column threshold to mydata 
mydata=merge(mydata,new_thr,by="seq_names",all=TRUE)




#verify monophyly of clusters

##create data.frame with leaves names and corresponding nodes in mytree$edge[,2] 
mytree$edge
tip_N=NULL
tip_edge=NULL
for (x in 1:length(mytree$tip.label))
{
	tip_edge=c(tip_edge,which.edge(mytree,mytree$tip.label[x]))
	#print(tip_edge)
	tip_N=c(tip_N,mytree$edge[tip_edge[x],2])
	#print(mytree$tip.label[x])
	#print(tip_nodes)
	}
tip_nodes=data.frame(mytree$tip.label,tip_N)
names(tip_nodes)=c("seq_names","tree_node")

##add column with tree nodes for sequences(tips)
mydata=merge(mydata,tip_nodes,by="seq_names",all=TRUE)
names(mydata)=c(names(mydata)[1:5],"thr_0.01","thr_0.02","thr_0.05","thr_0.07","thr_0.1",names(mydata)[11])
##problems with factros and numbers as characterc in the matrix mydata solved by passing to dataframe and changing to characters and integers
mydata.df=as.data.frame(mydata)
names(mydata.df)=c(names(mydata.df)[1:5],"thr_0.01","thr_0.02","thr_0.05","thr_0.07","thr_0.1",names(mydata.df)[11])
for (x in 1:5)
{
	mydata.df[,x]=as.vector(mydata.df[,x])
	}
for (x in 6:11)
{
	mydata.df[,x]=as.integer(mydata.df[,x])
	}
write.table(mydata.df,"CO1data_w_clusters.txt",sep="\t")




#renaming internal nodes labels as the nodes appear in the tree and giving node_names to be printed in the tree
mytree$node.label=((length(mytree$tip)+1):((length(mytree$tip)*2)-1))
aux=NULL
for (x in 1:length(mytree$node.label))
{
	aux[x]=paste("n",x,sep="")
	#print(aux)
	}
internal_nodes=data.frame(aux,mytree$node.label)
names(internal_nodes)=c("node_label","tree_node")
write.table(internal_nodes,"renaming_intNodes.txt",sep="\t")
#pdf('mytree_nodes.pdf')
#plot.phylo(mytree,show.tip.label=TRUE,show.node.label=TRUE,tip.color="purple",cex=0.5)
#dev.off()




##verifying the monophyly of clusters for each threshold
aux_max=NULL
for (x in 6:10)
	{
	aux_max=c(aux_max,max(mydata.df[,x]))
	}
	aux_max=as.integer(aux_max)
monophyl=matrix(NA,ncol=length(thr),nrow=max(aux_max), dimnames=list(c(1:max(aux_max)),thr))
for (x in 6:10)
	{
	u=unique(mydata.df[,x])   #vector with all the different clusters from columns 6 to 10 have the cluster results for the various thresholds
	for (y in u)
		{
		clade=mydata.df[mydata.df[,x]==y,1] #gets all the lines in column 1 that have in columnn x the value y, ie, all labels for cluster y
		monophyl[y,x-5]=is.monophyletic(mytree,clade)
		if (monophyl[y,x-5]=="FALSE") 
			{
			nonMono_clades[[paste(y,names(mydata.df)[x],sep="_")]]=clade     #vectors' names in nonMono_clades list correspond to the row.names in aux_ClusNonMono
			}
		}
	}
save(nonMono_clades,file="nonMono_clades.RData")   #to save the list as a file .RData that can be loaded alone
write.table(summary(nonMono_clades),"summary clades non_monophyl.txt")
colnames(monophyl)=c("thr_0.01","thr_0.02","thr_0.05","thr_0.07","thr_0.1")
write.table(monophyl,"monophyly_test.txt",sep="\t")




#create statistics of non-monophyletic groups per threshold
non_monofi=matrix(NA,5,2)
rownames(non_monofi)=c("thr_0.01","thr_0.02","thr_0.05","thr_0.07","thr_0.1")
colnames(non_monofi)=c("false","true")
for (x in 1:5)
{
non_monofi[x,]=table(monophyl[,x])
}
no_monofi=cbind(non_monofi,non_monofi[,1]+non_monofi[,2])
write.table(no_monofi,"stat_non_monophy.txt",sep="\t",col.names=c("non_monophy","monophy","Cl_total"))



#correct non-monophyletic groups

##for sequences to add in non monophy groups: get cluster number on column 6:10 (depending on the thr in nonMono_cluster name) of line with a given value in column4 (indiv_id)
seqs_forMono=read.table("seqs_for_monophy.txt", header=TRUE)
clust_old_nb=rep(NA,length(seqs_forMono[,1]))
seqs_forMono=cbind(seqs_forMono,clust_old_nb)
for(x in 1:length(seqs_forMono[,1]))
{
y=as.character(seqs_forMono[x,2])
seqs_forMono[x,4]=mydata.df[mydata.df[,4]==seqs_forMono[x,1],names(mydata.df)==y]
}
write.table(seqs_forMono,"seqs to move to nonMonophy groups.txt")

##break seqs_forMono table in subtables per threshold (analysis will be easier when comparing with mydata.df)
thr001=seqs_forMono[seqs_forMono[,2]=="thr_0.01",]
thr002=seqs_forMono[seqs_forMono[,2]=="thr_0.02",]
thr005=seqs_forMono[seqs_forMono[,2]=="thr_0.05",]
thr007=seqs_forMono[seqs_forMono[,2]=="thr_0.07",]
thr01=seqs_forMono[seqs_forMono[,2]=="thr_0.1",]

##convert cluster numbers in thr00x into factors and count frequency of each cluster 
old_clu001=as.data.frame(table(as.factor(thr001[,4])))
old_clu002=as.data.frame(table(as.factor(thr002[,4])))
old_clu005=as.data.frame(table(as.factor(thr005[,4])))
old_clu007=as.data.frame(table(as.factor(thr007[,4])))
old_clu01=as.data.frame(table(as.factor(thr01[,4])))
names(old_clu001)[1]="old_cluster"
names(old_clu002)[1]="old_cluster"
names(old_clu005)[1]="old_cluster"
names(old_clu007)[1]="old_cluster"
names(old_clu01)[1]="old_cluster"

##count the frequency of each cluster per threshold in the mydata.df table
mydata.df001=as.data.frame(table(mydata[,6]))
mydata.df002=as.data.frame(table(mydata[,7]))
mydata.df005=as.data.frame(table(mydata[,8]))
mydata.df007=as.data.frame(table(mydata[,9]))
mydata.df01=as.data.frame(table(mydata[,10]))
mydataClu001=NULL
mydataClu002=NULL
mydataClu005=NULL
mydataClu007=NULL
mydataClu01=NULL

for (x in 1:length(old_clu001[,1]))
{
mydataClu001=c(mydataClu001,mydata.df001[as.character(mydata.df001[,1])==as.character(old_clu001[x,1]),2])
}
for (x in 1:length(old_clu002[,1]))
{
mydataClu002=c(mydataClu002,mydata.df002[as.character(mydata.df002[,1])==as.character(old_clu002[x,1]),2])
}
for (x in 1:length(old_clu005[,1]))
{
mydataClu005=c(mydataClu005,mydata.df005[as.character(mydata.df005[,1])==as.character(old_clu005[x,1]),2])
}
for (x in 1:length(old_clu007[,1]))
{
mydataClu007=c(mydataClu007,mydata.df007[as.character(mydata.df007[,1])==as.character(old_clu007[x,1]),2])
}
for (x in 1:length(old_clu01[,1]))
{
mydataClu01=c(mydataClu01,mydata.df01[as.character(mydata.df01[,1])==as.character(old_clu01[x,1]),2])
}
old_clu001=cbind(old_clu001,mydataClu001)
old_clu002=cbind(old_clu002,mydataClu002)
old_clu005=cbind(old_clu005,mydataClu005)
old_clu007=cbind(old_clu007,mydataClu007)
old_clu01=cbind(old_clu01,mydataClu01)
write.table(old_clu001,"stat_seqs_to_cluster_into monophyl0.01.txt")
write.table(old_clu002,"stat_seqs_to_cluster_into monophyl0.02.txt")
write.table(old_clu005,"stat_seqs_to_cluster_into monophyl0.05.txt")
write.table(old_clu007,"stat_seqs_to_cluster_into monophyl0.07.txt")
write.table(old_clu01,"stat_seqs_to_cluster_into monophyl0.1.txt")####compare columns in old_clu00xto check if sequences to be moved will make an old cluster become part of a new one 

##changing old_cluster to new cluster in mydata.df
mydata.dfmonofi=mydata.df
for(x in 1:length(thr001[,4]))
{mydata.dfmonofi[mydata.dfmonofi[,6]==thr001[x,4],6]=thr001[x,3]}
for(x in 1:length(thr002[,4]))
{mydata.dfmonofi[mydata.dfmonofi[,7]==thr002[x,4],7]=thr002[x,3]}
for(x in 1:length(thr005[,4]))
{mydata.dfmonofi[mydata.dfmonofi[,8]==thr005[x,4],8]=thr005[x,3]}
for(x in 1:length(thr007[,4]))
{mydata.dfmonofi[mydata.dfmonofi[,9]==thr007[x,4],9]=thr007[x,3]}
for(x in 1:length(thr01[,4]))
{mydata.dfmonofi[mydata.dfmonofi[,10]==thr01[x,4],10]=thr01[x,3]}
write.table(mydata.dfmonofi,"CO1data_w_clusters_monophy.txt")

##verifying the monophyly in mydata.dfmonofi clusters
aux_maxfi=NULL
nonMono_cladesfi=list()
for (x in 6:10)
	{
	aux_maxfi=c(aux_maxfi,max(as.numeric(mydata.dfmonofi[,x])))
	}
monophylfi=matrix(NA,ncol=length(thr),nrow=max(aux_maxfi), dimnames=list(c(1:max(aux_maxfi)),thr))
for (x in 6:10)
	{
	u=unique(mydata.dfmonofi[,x])
	for (y in 1:length(u))
		{
		clade=mydata.dfmonofi[mydata.dfmonofi[,x]==u[y],1] 
		#print(clade)
		monophylfi[y,x-5]=is.monophyletic(mytree,clade)
		if (monophylfi[y,x-5]=="FALSE") 
			{
			nonMono_cladesfi[[paste(u[y],names(mydata.dfmonofi)[x],sep="_")]]=clade    
			}
		}
	}
print(nonMono_cladesfi)
save(nonMono_cladesfi,file="retest_nonMono_cladesfi.RData")  
#write.table(summary(nonMono_clades),"summary clades non_monophyl.txt")
#print(monophyl)
colnames(monophylfi)=c("thr_0.01","thr_0.02","thr_0.05","thr_0.07","thr_0.1")
write.table(monophylfi,"monophyly_retest.txt",sep="\t")

################
#checking which clusters have individuals from 3 or 4 islands
mydata.dfmonofi[,2]=as.factor(mydata.dfmonofi[,2])
for (y in 6:10)
{mydata.dfmonofi[,y]=as.numeric(mydata.dfmonofi[,y])}
for (x in 6:10)
{
	isl3to4=data.frame()
	clust=unique(mydata.dfmonofi[,x])
	for (w in clust)
	{
		subtable=mydata.dfmonofi[mydata.dfmonofi[,x]==w,c(1:5,x)]		#print(subtable)
		if(length(unique(subtable$island))>=3)
		{
			isl3to4=rbind(isl3to4,subtable)
		}
	}
	
	write.table(isl3to4,paste("island3to4",names(mydata.dfmonofi[x]),sep="_"))
}

island_001=read.table("island3to4_thr_0.01")
island_002=read.table("island3to4_thr_0.02")
island_002=island_002[,-6]
island_005=read.table("island3to4_thr_0.05")
island_002=island_002[,-(6:7)]
island_007=read.table("island3to4_thr_0.07")
island_002=island_002[,-(6:8)]
island_01=read.table("island3to4_thr_0.1")
island_002=island_002[,-(6:9)]





##some statistics on groups with >=3 islands per threshold
vector001=as.factor(island_001$thr_0.01)
table(vector001)
count_gr=length(levels(vector001))
vector002=as.factor(island_002$thr_0.02)
table(vector002)
count_gr=cbind(count_gr,length(levels(vector002)))
vector005=as.factor(island_005$thr_0.05)
table(vector005)
count_gr=cbind(count_gr,length(levels(vector005)))
vector007=as.factor(island_007$thr_0.07)
table(vector007)
count_gr=cbind(count_gr,length(levels(vector007)))
vector01=as.factor(island_01$thr_0.1)
table(vector01)
count_gr=cbind(count_gr,length(levels(vector01)))
colnames(count_gr)=c("NbCluster001","NbCluster002","NbCluster005","NbCluster007","NbCluster01")
write.table(count_gr,"number of clusters_thr for subtrees.txt")
xp="NA"
xp001=rep(xp,24)
xp002=rep(xp,20)
xp005=rep(xp,13)
xp007=rep(xp,11)
vector001=c(levels(vector001),xp001)
vector002=c(levels(vector002),xp002)
vector005=c(levels(vector005),xp005)
vector007=c(levels(vector007),xp007)
vector01=levels(vector01)
gr_sub3=cbind(vector001,vector002,vector005,vector005,vector007,vector01)
write.table(gr_sub3,"groups per thr for subtrees.txt")


######implement function(threshold) to accept the vector of 5 thresholds:0.01,0.02,0.05,0.07 or 0.1



##assigning node numbers as labels for internal nodes
nodeLabels(my3_phylo4)=paste("N", nodeId(my3_phylo4, "internal"),sep"")

save.image("CO1_production of groups and subtrees.RData")